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A recently developed lattice method for large numbers of strongly interacting nonrelativistic 
fermions exhibits a heavy tail in the distributions of correlators for large Euclidean time T and 
large number of fermions N, which only allows the measurement of ground state energies for a 
limited number of fermions using standard techniques. In such cases, it is suggested that mea- 
suring the log of the correlator is more efficient, and a cumulant expansion of this quantity can 
be exactly related to the correlation function. The cumulant expansion technique allows us to 
determine the ground state energies of up to 66 unpolarized unitary fermions on lattices as large 
as 72 x 14 3 , and up to 70 unpolarized unitary fermions trapped in a harmonic potential on lattices 
as large as 72 x 64 3 . We have also improved our lattice action with a Galilean invariant form for 
the four-fermion interaction, which results in predictive volume scaling of the lowest energy of 
three fermions in a periodic box and in good agreement of our results for N <6 trapped unitary 
fermions with those from other benchmark calculations. 
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1. Introduction 



One of the most challenging problems in physics is the quantitative understanding of a quan- 
tum system of many strongly interacting particles for which numerical simulations have played an 
important role. Among the variety of known systems, a dilute Fermi gas near unitarity has been 
noted for its pure form, which can be readily studied from theory, and for its realization by using 
ultra-cold atomic experiments (for a recent review, see [1, 2] ). Besides the intrinsic physical inter- 
est, such as universal behavior, the unitary Fermi gas might also be considered as an ideal starting 
point to develop numerical techniques which can be applied to low energy nuclear physics, and to 
attack the noise problem which typically appears in numerical simulations of many particles. 

In 2010 we presented a highly improved lattice method for non-relativistic fermions with four- 
fermion contact interactions [3] and its applications for trapped and untrapped unitary fermions 
[4, 5], where the largest numbers of fermions were restricted to N = 20 and 38, respectively, due 
to a statistical overlap problem. Since then, we have devised a tuning technique with a Galilean 
invariant form for the four-fermion interaction, which keeps the system at unitary up to Galilean 
boosts. We have also implemented an external harmonic potential in a more sophisticated way. 
Finally, we have developed a cumulant expansion technique to extract the ground state energies 
from data exhibiting a distribution overlap problem [6, 7]. As a result, we are able to extend our 
calculation of ground state energies of unpolarized fermions up to ,/V = 70 and up to /V = 66 unitary 
fermions with and without a harmonic trap, respectively. For the untrapped case, this extension 
leads us to the regime, N > 38, where the ground state energy E in units of that for non-interacting 
particles 1 is constant; this implies that we are near the thermodynamic limit and can determine 
the universal dimensionless parameter t, = E/E^°\ called the Bertsch parameter. On the other 
hand, for the trapped case we find that N = 70 is insufficient to reach the thermodynamic limit. 



2. Lattice construction for fermions at unitarity 

We consider simulations of N non-relativistic two-component fermions \j/ = on a 

T x L 3 Euclidean lattice, where periodic boundary conditions on the spatial directions and an open 
boundary condition on the temporal direction are imposed. The four-fermion contact interaction 
between different species of fermions is generated by Z 2 auxiliary fields which live on the time- 
like links of the lattice. A consequence of this set up is that the fermion-determinant does not 
depend on the <p field and thus the quenched simulation is exact [3] . 

The discretization errors for the energy of a single fermion are eliminated through the use of 
a perfect dispersion relation, while the unitarity limit may be achieved by tuning the interaction 
strength C(q) to reproduce scattering data, pcot^o = 0, using the Liischer formula [3, 7]. Here 
the operator q for |q| < A corresponds to the momentum transfer between incoming and outgoing 
fermions, which makes the interaction Galilean invariant. The periodic function C(q) is defined by 

C ( t ») = 17 I c 2„0 2 „(q), (2.1) 



'The energies of non-interacting untrapped and trapped fermions are given by ^^trapped = (3A0^ 5//3 ^ 4/ ^/lOML 



and ^tfjjppg^ = (3iV) 4 / 3 fi)/4, respectively. Here, CO is an oscillator frequency. 
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Figure 1: ln-ln plot of p cot So (dots) along with 
expected 77 scaling (dashed lines) for which No = 
1, 2, 3 and 4 coefficients are tuned by exactly match- 
ing the first No energy eigenvalues for two particles 
in a box. Data is from an L = 32 and M = 5 lattice. 



Figure 2: L dependence of the 9th energy eigenval- 
ues for two particles in a box. The dashed lines are 
fit results of the data (dots) using Eq. 2.3. 



where the basis functions are 



2n (q) 



Ml x < 



. e -1 2 /M 



|q| < A 

|q| > A' 



(2.2) 



for q within the first Brillouin zone and periodic from one Brillouin zone to the next. Although the 
parameter Mq may in general be different from the fermion mass M, we use Mq = M in this work. 

The coefficients C2 n are numerically determined by matching the lowest No energy eigenvalues 
of the two-body transfer matrix in our lattice theory with the lowest No Luscher energy eigenvalues 
in the corresponding continuum theory. In Fig. 1 we plot the /?cot5o (dots), computed using the 
exact lattice energy eigenvalues and Luscher's formula, which results from tuning the first No terms 
in the effective range expansion to zero. These results show an expected 77 scaling represented by 
the dashed lines in the figure. A systematic improvement is also seen for the energies of excited 
states above the lowest No tuned states, where the correction may be predicted by [7], 



\\ocL 



2-2n 



(2.3) 



Here 7]^ are the kth eigenvalues of the two-body transfer matrix with No terms tuned, while 77^ are 
the kth solutions of Luscher's formula in the unitary limit. As an example, we plot L(j\^/i]t — 1) 



with respect to L (dots) in Fig. 2 and find good agreement with Eq. 2.3. 

In simulations of trapped unitary fermions, the external harmonic potential, U 
spring constant jc, has been implemented in the transfer matrix as 



with 



tapped = e^^e^iX - b % V)e- b * u e- h ^l m 



-b z {je+u)+o{b>) 



(2.4) 



where V represents the interaction and Ji? + U is the target Hamiltonian for trapped unitary fermions. 
As seen in this formula, temporal discrerization errors appear at 0{b 2 z ). 2 



2 Throughout this paper, b s and b x represent the spatial and temporal lattice spacings, respectively. 
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3. Measurement, overlap problem, and cumulant expansion method 

The jty-body correlators are constructed by taking determinants of the Slater matrix, where each 
element is obtained by evolving an initial state (source) at t = with a single particle propagator 
and projecting onto a final state (sink) at z = T. The choice of sources and sinks for untrapped and 
trapped unitary fermions, which gives superior wave function overlap by imposing two-particle 
correlations at the sinks, has been described in [4] and [5], respectively. 

One of the greatest difficulties in extending our study to large N was the apparent upward 
drift of the effective mass plot at large Euclidean time. By investigating the distribution of A^-body 
correlators, we found that this difficulty arises from a heavily long-tailed distribution, requiring 
an exponentially large number of samples before the central limit theorem becomes applicable. 
In other words, the path-integral probability measure has small overlap with the dominant part 
of the operator being estimated. On the other hand, the distribution of the log of correlators is 
nearly Gaussian, implying that the standard estimation of the lower moments should succeed with 
moderately sized ensembles. A general relation between the expectation value of & and that of 
ln^ 7 is given by [6] 

Mn*)) = t*, (3-D 

n=\ "• 

where K k is the kth cumulant of ln'rf . The generalized effective mass and the ground state energy 
associated with each partial sum in Eq. 3.1 may be written as [7] 

1 Nk 1 

m N e ff^) = -^T,-[K^ + Ar)-K n (z)} and E Nk = lim mj f (t). (3.2) 

Since the statistical uncertainties typically increase as A^ increases, one may determine the ideal 
value N£ for which the statistical uncertainties and truncation errors are comparable. 

4. Ground state energy of unitary fermions 

To extract the energies of the system, we perform standard bootstrap resamplings and corre- 
lated % 2 fits to the plateau region of the generalized effective mass. Fitting systematic errors are 
obtained by varing the end points of the fitting interval. The quoted errors represent the combina- 
tion of the statistical and fitting systematic errors added in quadrature. 

4.1 Untrapped unitary fermions 

As a nontrivial test of our lattice method, we have computed the lowest energy of three unitary 
fermions in a zero total momentum state with high precision. We performed the calculation for 
lattice sizes L = 8, 10, 12, 14, 16, tuning the coefficients of four C?2n operators for the L = 8 lattice, 
and five for the other lattices. Since the single particle and two particle s-wave sectors have been 
highly tuned, the leading L dependence will be L~ 3 due to the untuned two-derivative two-body 
p-wave operator. In Fig. 3 we plot the lowest energies versus L~ 3 . We perform two parameter fits 
of the L > 10 data to c\ + C2/L 3 and find an infinite volume energy of 0.3735^qqqq7 in units of 
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Figure 3: Energy of three untrapped unitary 
fermions in a zero total momentum eigenstate 
plotted versus (b s /L) 3 for L/b s = 8, 10, 12, 14, 16. 
The red band represents the uncertainty in two- 
parameter fits of the L/b s > 10 data to the function 
c\ +C2/L 3 . while the black line is the fit to the cen- 
tral values. The dashed line is the result from Ref. 
[8]. 
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Figure 4: (Preliminary) Plot of ground state ener- 
gies versus numbers of untrapped unitary fermions 
for L = 10, 12 and 14. For L = 12 we considered 
three different sinks, j3 = 0.5,0.75, 1.0, to take into 
account any source dependence in our extraction of 
the energies. The arrows indicate the positions for 
which noninteracting fermions have fully occupied 
a given momentum shell. 



the energy of three noninteracting fermions, which agrees with the high precision calculation by 
Pricoupenko and Castin [8] within our ~ 0.3% uncertainty. 

The preliminary results for the ground state energies of up to 66 unpolarized unitary fermions 
in a periodic box are shown in Fig. 4. With the given statistics, we do not resolve any shell structure 
in the energy for N > 38 and believe that the system is close to the thermodynamic limit. By 
averaging the results obtained from a constant fit of B, (N) for TV > 38 using five different ensembles, 
we find a preliminary result for the Bertsch parameter, E, = 0.399 ±0.002. The unitary Fermi gas 
has been extensively studied using the Quantum Monte Carlo (QMC) technique and the most recent 
calculation reported is E, < 0.383(3) [9, 10]. Recent experimental results for this parameter are 
§ =0.39(2) [11] and 0.41(1) [12]. 

4.2 Trapped unitary fermions 

In numerical simulations of the trapped unitary Fermi gas, we have additional time scale 1 / (O 
and length scale Lq = (Moo) -1 / 2 which should be chosen so that b s <C Lo <C L and b x <C l/co 
in order to minimize the discretization and finite volume errors. To balance the need for small 
temporal discretization errors with the computational cost associated with the number of time steps 
required to reach the ground state, we have chosen (ob x = 0.005. An ideal Lo has been determined 
by scanning the parameter space of Lo/b s and L/Lo and finding the region where both finite volume 
and spatial discretization errors are small. Fig. 5 presents our findings for the ground state energies 
of N < 6 fermions, with Lo/b s ranging from 3 to 8 and fixed L/b„ = 48. For Lo/b s < 7 we find 
that the systematic errors increase as Lo decreases, which indicates that the discretization error is 
not negligible. For Lo/b s > 7 we have also performed simulations with L/b s = 64, which showed 
no volume dependence; and thus we conclude that both types of systematic errors are negligible. 
For 6 < N < 70, based on the results of the Lo scan for TV < 6, we have chosen to perform the 
calculation for three volumes (L = 48,54,64) and at two values of the trap size (Lq = 7.5,8). We 
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Table 1: Results for^trapped/' 3 for /V < 6, including combined statistical and fitting systematic errors (first 
row). For comparison we give the exact N = 3 result [13] and results of Ref. [14] (second and third rows). 
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Figure 5: Ground-state energies (in units of co) as 
a function of Lo/b s at fixed L/b s = 48 for various 
values of N. Dashed lines are results from Ref. [14]. 



Figure 6: Ground-state energies of Af-trapped uni- 
tary fermions. 



find non-negligible volume dependence and perform an infinite volume extrapolation for each Lq. 
All of these considerations are included along with statistical, fitting systematic, and truncation of 
the cumulant expansion errors in the final quoted errors of ground state energies. 

We first benchmark our method for up to N = 6 against high-precision solutions to the many- 
body Schrodinger equation [14], achieving agreement at 1% as shown in Tab. 1. In Fig. 6 we 
plot the results of the ground state energies for N < 70 along with the results from two fixed- 
node calculations for comparison: a Green's function Monte Carlo (GFMC) approach [15] and a 
diffusion Monte Carlo (FN-DMC) approach [16], which provide upper bounds on the ground state 
energies. We find that our energies are consistently lower than those obtained using both of these 
methods. However, our results show clear shell structure which indicates that^V ~ 70 is insufficient 
to reach the thermodynamic limit. 



5. Conclusion 

A new statistical technique, called the cumulant expansion method, allows us to extend our 
studies up to ./V ~ 70 unitary fermions both with and without an external harmonic trap. We re- 
port a preliminary value of ^ = 0.399 ±0.002 from the calculation of the ground state energy for 
N >3S untrapped unitary fermions. For N <10 trapped fermions, we find that the shell structure 
is pronounced, which implies that we have not reached the thermodynamic limit, but our values 
of the ground state energies are consistently lower than those obtained from variational QMC cal- 
culations. In order to check the validity of our lattice method, we have performed high-precision 
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calculations for N = 3 untrapped fermions and N <6 trapped fermions, and found good agreement 
with the results in [8] and [13, 14], respectively. 

In the near future, we will present results for the pairing gap by calculating the ground state 
energies of the slightly unpolarized unitary Fermi gas (A^ = + 1) and the integrated contact 
density by studying the dependence of the ground state energies on the s-wave scattering length near 
the unitary regime. Our lattice method combined with the new statistical technique is applicable 
for a wide range of nonrelativistic many-body systems. 
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